Load in Data

load("./data/elk_processed.rda")
library(sf)

library(dplyr)
elk_res <- elk_processed |>
  filter(id %in% c("YL58","YL64","YL80", "YL91", "YL94"))
elk_res_sf <- elk_res |>
  st_as_sf(coords = c("lon","lat"), crs=4326) |>
  st_transform(32611)
trails_and_roads <- st_read("./Hebblewhitematerials/data/humanacess.shp") |>
  st_transform(32611)
## Reading layer `humanacess' from data source 
##   `C:\Users\nicol\Documents\BrazilMove2024\Hebblewhitematerials\data\humanacess.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 24134 features and 39 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 513240.9 ymin: 5665028 xmax: 650430.4 ymax: 5760389
## Projected CRS: NAD83 / UTM zone 11N
str(trails_and_roads)
## Classes 'sf' and 'data.frame':   24134 obs. of  40 variables:
##  $ IGDS_LAYER: chr  NA NA NA NA ...
##  $ IGDS_LEVEL: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ IGDS_COLOR: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ NAME      : chr  NA NA "Spreading/ Porcupine conn" "Porcupine Cr. Tr." ...
##  $ CODE      : chr  "DC31700100" "DC31700100" "DC31700100" "DC31700100" ...
##  $ TYPE      : chr  "TRAIL" "Trail" "Trail" "Trail" ...
##  $ CLASS     : chr  NA NA NA NA ...
##  $ SOURCE    : chr  "BAG" "OldHUM" "OldHUM" "OldHUM" ...
##  $ ACCURACY  : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ STATUS    : chr  NA NA NA NA ...
##  $ NOTES     : chr  NA NA NA NA ...
##  $ MOTOR     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ OHV       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ FOOT      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ CYCLE     : int  0 0 0 0 0 0 0 0 1 1 ...
##  $ HORSE     : int  0 0 0 0 0 0 1 1 1 1 ...
##  $ SKI       : int  0 0 0 0 0 0 0 0 1 1 ...
##  $ SNOWM     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ COMMERCIAL: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ OTHER     : chr  NA NA NA NA ...
##  $ SUM_CLASS : chr  "MEDIUM" "MEDIUM" "MEDIUM" "MEDIUM" ...
##  $ SUM_USE   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ESTIMATE  : chr  "BAG" "BAG" "BAG" "BAG" ...
##  $ CONFIDENCE: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ SUM_YEAR  : chr  "2001" "2002" "2002" "2002" ...
##  $ WIN_CLASS : chr  "0" NA NA NA ...
##  $ WIN_USE   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ ESTIMATE0 : chr  "0" "0" "0" "0" ...
##  $ CONFIDEN0 : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ WIN_YEAR  : chr  "0" "0" "0" "0" ...
##  $ MAP_SHEET : chr  NA NA NA NA ...
##  $ YEAR_MAP  : chr  NA "2001" "2001" "2001" ...
##  $ YEAR_UPD  : chr  "2002" "2002" "2002" "2002" ...
##  $ DISPOSITIO: chr  NA NA NA NA ...
##  $ SOURCETHM : chr  NA NA NA NA ...
##  $ MOTO_YEAR : chr  NA NA NA NA ...
##  $ YEAR      : chr  NA NA NA NA ...
##  $ COMMENT   : chr  NA NA NA NA ...
##  $ LENGTH    : num  1013 1469 4149 4079 3988 ...
##  $ geometry  :sfc_MULTILINESTRING of length 24134; first list element: List of 1
##   ..$ : num [1:6, 1:2] 543862 544276 544400 544438 544469 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTILINESTRING" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:39] "IGDS_LAYER" "IGDS_LEVEL" "IGDS_COLOR" "NAME" ...
st_crs(trails_and_roads)==st_crs(elk_res_sf)
## [1] TRUE
plot(st_geometry(trails_and_roads))
plot(elk_res_sf$geometry, add=TRUE, col="green")

library(raster)
landcover <- raster("./Hebblewhitematerials/data/landcover.tif")
landcover
## class      : RasterLayer 
## dimensions : 3811, 2652, 10106772  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 540827.8, 620387.8, 5642002, 5756332  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs 
## source     : landcover.tif 
## names      : TRUE. 
## values     : 0, 16  (min, max)
plot(landcover)

landcover <- projectRaster(landcover, crs = crs(elk_res_sf))

0 - NA 1 - Open Conifer Forest 2 - Moderate Conifer Forest 3 - Closed Conifer Forest 4 - Deciduous Forest 5 - Mixed Forest 6 - Regeneration 7 - Herbaceous 8 - Shrub 9 - Water 10 - Rock-Ice 11 - Cloud 12 - Burn-Forest 13 - Burn-Grassland 14 - Burn-Shrub 15 - Alpine Herb 16 - Alpine Shrub

sort(unique(round(landcover@data@values)))
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
landcover_rc <- reclassify(landcover, c(-Inf, 0, 0,
                                       0, 1, 1,
                                       1, 2, 1,
                                       2, 3, 1,
                                       3, 4, 1,
                                       4, 5, 1,
                                       5, 6, 2,
                                       6, 7, 3,
                                       7, 8, 3,
                                       8, 9, 4,
                                       9, 10, 5,
                                       10, 11, 6,
                                       11, 12, 7,
                                       12, 13, 7,
                                       13, 14, 7,
                                       14, 15, 8,
                                       15, Inf, 8))

0 - NA 1 - Forest (Conifer, Deciduous, Mixed) 2 - Regeneration 3 - Herbaceous/Shrub 4 - Water 5 - Rock-Ice 6 - Cloud 7 - Burned (Forest, Grassland, Shrub) 8 - Alpine (Herb, Shrub)

sort(unique(round(landcover_rc@data@values)))
## [1] 0 1 2 3 4 5 6 7 8
plot(landcover_rc)
plot(st_geometry(elk_res_sf), add=TRUE)

Resource Selection Function

Determine “Available” Locations

library(sp)

library(adehabitatHR)
elk_mcp <- mcp(as_Spatial(elk_res_sf), 100)
determineAvailableLocations <- function(used_data, mcp){
  
  random_n <- nrow(used_data)+1000
    
  sample <- spsample(mcp, random_n, "random") |> st_as_sf()
  
  sample$id <- used_data$id[1]
  
  sample$Used <- FALSE
  used_data$Used <- TRUE
  
  sample <- sample |> dplyr::select(id, Used, geometry)
  used_data <- used_data |> dplyr::select(id, Used, geometry)
  
  combine_data <- rbind(used_data, sample)
  
  return(combine_data)
}
elk_id <- split(elk_res_sf, elk_res_sf$id)
elk_id_combine <- lapply(elk_id,
                         determineAvailableLocations,
                         elk_mcp)

elk_combine <- do.call("rbind", elk_id_combine)

Create Covariates

elk_mcp_sf <- elk_mcp |>
  st_as_sf()
trails_and_roads2 <- st_crop(trails_and_roads, elk_mcp_sf)

do elk come within 250 m of trails?

trails_and_roads_buff <- st_buffer(trails_and_roads2, 250) 
library(mapview)
mapview(trails_and_roads2)+
  mapview(trails_and_roads_buff, col.regions="red")+
  mapview(subset(elk_combine, Used==TRUE))
trails_and_roads_buff$Near_Trail <- TRUE

trails_and_roads_buff <- trails_and_roads_buff[, c(40:41)]

str(trails_and_roads_buff)
## Classes 'sf' and 'data.frame':   513 obs. of  2 variables:
##  $ geometry  :sfc_POLYGON of length 513; first list element: List of 1
##   ..$ : num [1:232, 1:2] 590607 590566 590561 590528 590517 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
##  $ Near_Trail: logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
##   ..- attr(*, "names")= chr "Near_Trail"
elk_combine2 <- st_join(elk_combine, trails_and_roads_buff, left=TRUE)

elk_combine2[which(is.na(elk_combine2$Near_Trail)),]$Near_Trail <- FALSE
mapview(trails_and_roads2)+
  mapview(subset(elk_combine2, Used==TRUE), zcol="Near_Trail")

Extract Resource (Raster) Layer to Points

elk_combine_sp <- as(elk_combine2, Class="Spatial")
landcover_elk <- extract(landcover_rc, elk_combine_sp)
elk_combine2$Landcover <- landcover_elk

Prepare Data for Modeling

0 - NA 1 - Forest (Conifer, Deciduous, Mixed) 2 - Regeneration 3 - Herbaceous/Shrub 4 - Water 5 - Rock-Ice 6 - Cloud 7 - Burned (Forest, Grassland, Shrub) 8 - Alpine (Herb, Shrub)

elk_df <- data.frame(elk_combine2)

elk_df$Landcover <- as.factor(elk_df$Landcover)

elk_df2 <- elk_df |> dplyr::mutate(No_Data = grepl("0", Landcover),
                         Forest = grepl("1", Landcover),
                         Regeneration = grepl("2", Landcover),
                         Herbaceous_Shrub = grepl("3", Landcover),
                         Water = grepl("4", Landcover),
                         Rock_Ice = grepl("5", Landcover),
                         Cloud = grepl("6", Landcover),
                         Burned = grepl("7", Landcover),
                         Alpine = grepl("8", Landcover))

Fit RSF

library(ggplot2)
landcover_count <- count(elk_df2, Landcover, Used)

ggplot() +
  geom_bar(data = landcover_count, aes(x = Landcover, y = n, fill = Used), position="dodge", stat="identity") +
  scale_x_discrete(labels=c("Forest","Regeneration","Herbaceous_Shrub","Water", "Rock_Ice","Cloud","Burned","Alpine")) +
  theme_classic() +
  labs(x = "Landscape Class", y = "Count")

Binomial GAM

elk_glm <- glm(Used ~ Near_Trail + Forest + Herbaceous_Shrub + Water + Burned + Alpine, data = elk_df2, family = "binomial")

summary(elk_glm)
## 
## Call:
## glm(formula = Used ~ Near_Trail + Forest + Herbaceous_Shrub + 
##     Water + Burned + Alpine, family = "binomial", data = elk_df2)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.66391    0.02636 -63.121  < 2e-16 ***
## Near_TrailTRUE        1.09660    0.01823  60.154  < 2e-16 ***
## ForestTRUE            0.13086    0.02818   4.643 3.43e-06 ***
## Herbaceous_ShrubTRUE  2.49255    0.02831  88.058  < 2e-16 ***
## WaterTRUE             1.61280    0.06488  24.858  < 2e-16 ***
## BurnedTRUE           -0.92736    0.05506 -16.842  < 2e-16 ***
## AlpineTRUE           -0.25616    0.05918  -4.329 1.50e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 131839  on 96462  degrees of freedom
## Residual deviance:  93102  on 96456  degrees of freedom
## AIC: 93116
## 
## Number of Fisher Scoring iterations: 4
elk_glm2 <- glm(Used ~ Near_Trail + Herbaceous_Shrub + Water + Burned + Alpine, data = elk_df2, family = "binomial")

summary(elk_glm2)
## 
## Call:
## glm(formula = Used ~ Near_Trail + Herbaceous_Shrub + Water + 
##     Burned + Alpine, family = "binomial", data = elk_df2)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.56973    0.01655 -94.870  < 2e-16 ***
## Near_TrailTRUE        1.10958    0.01801  61.608  < 2e-16 ***
## Herbaceous_ShrubTRUE  2.38852    0.01714 139.357  < 2e-16 ***
## WaterTRUE             1.51009    0.06099  24.761  < 2e-16 ***
## BurnedTRUE           -1.03080    0.05031 -20.487  < 2e-16 ***
## AlpineTRUE           -0.35774    0.05496  -6.509 7.54e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 131839  on 96462  degrees of freedom
## Residual deviance:  93124  on 96457  degrees of freedom
## AIC: 93136
## 
## Number of Fisher Scoring iterations: 4

no difference because there are no alpine values

anova(elk_glm, elk_glm2, test="Chi")
## Analysis of Deviance Table
## 
## Model 1: Used ~ Near_Trail + Forest + Herbaceous_Shrub + Water + Burned + 
##     Alpine
## Model 2: Used ~ Near_Trail + Herbaceous_Shrub + Water + Burned + Alpine
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1     96456      93102                          
## 2     96457      93124 -1  -21.788 3.045e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The AIC (Akaike information criterion) is a measure of fit that penalizes for the number of parameters

want a lower AIC

AIC(elk_glm, elk_glm2)
##          df      AIC
## elk_glm   7 93116.38
## elk_glm2  6 93136.17
model_summary <- data.frame(summary(elk_glm)$coefficients)

glm_results <- model_summary |> 
          dplyr::mutate(low = Estimate - 2*Std..Error, high = Estimate + 2*Std..Error) 

glm_results$variable <- row.names(glm_results)

ggplot(glm_results, aes(x = Estimate, variable)) + 
  geom_errorbarh(aes(xmin = low, xmax = high)) + 
  geom_point() +
  theme_classic() +
  ylab("Landscape Class")+
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red", size=1.5) +
  scale_y_discrete(labels=c("Intercept","Alpine","Burned","Forest","Herbaceous_Shrub","Near_Trail","Water"))

Binomial GAMM

library(glmmTMB)
elk_glmm <- glmmTMB(Used ~ Near_Trail + Herbaceous_Shrub + Water + Burned + Alpine + (1|id), data = elk_df2, family = binomial)

summary(elk_glmm)
##  Family: binomial  ( logit )
## Formula:          
## Used ~ Near_Trail + Herbaceous_Shrub + Water + Burned + Alpine +      (1 | id)
## Data: elk_df2
## 
##      AIC      BIC   logLik deviance df.resid 
##  92986.6  93052.9 -46486.3  92972.6    96456 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.01207  0.1099  
## Number of obs: 96463, groups:  id, 5
## 
## Conditional model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.62547    0.05208  -31.21  < 2e-16 ***
## Near_TrailTRUE        1.10901    0.01805   61.45  < 2e-16 ***
## Herbaceous_ShrubTRUE  2.39455    0.01719  139.28  < 2e-16 ***
## WaterTRUE             1.51199    0.06102   24.78  < 2e-16 ***
## BurnedTRUE           -1.01921    0.05037  -20.24  < 2e-16 ***
## AlpineTRUE           -0.34604    0.05505   -6.29 3.26e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model_summary <- data.frame(summary(elk_glmm)$coefficients$cond)

glmm_results <- model_summary |> 
          dplyr::mutate(low = Estimate - 2*Std..Error, high = Estimate + 2*Std..Error) 

glmm_results$variable <- row.names(glmm_results)

ggplot(glmm_results, aes(x = Estimate, variable)) + 
  geom_errorbarh(aes(xmin = low, xmax = high)) + 
  geom_point() +
  theme_classic() +
  ylab("Landscape Class")+
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red", size=1.5) +
  scale_y_discrete(labels=c("Intercept","Alpine","Burned","Forest","Herbaceous_Shrub","Near_Trail","Water"))